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a^ 

A review is given of recent results about the computation of irrotational Darwin-Riemann 
configurations in general relativity. Such configurations are expected to represent fairly well 
the late stages of inspiralling binary neutron stars. 



§1. Introduction 



The classical Darwin ellipsoids are equilibrium figures of incompressible fluid 
\ bodies in a synchronous binary system*). Synchronous means that each body is 

■ spinning with respect to some inertial frame at the same angular velocity as the 

, orbital angular velocity. In this manner it always presents the same face to its com- 

panion. The vorticity of the fluid with respect to some inertial frame is then equal 
to twice the orbital angular velocity. Darwin-Riemann configurations ^-^ are general- 
Q> . izations of Darwin ellipsoids to arbitrary vorticity (i.e. non-synchronous spins). As 

Q \ detailed below, a subset of Darwin-Riemann ellipsoids is of particular importance for 

O^' the late stages of inspiralling binary neutron stars, which are expected to be one of 

the strongest sources of gravitational radiation for the interferometric detectors cur- 
rently under construction (GEO600, LIGO, TAMA and VIRGO). This subset is the 
irrotational Darwin-Riemann configurations, i.e. configurations for which the fluid 
^ , vorticity with respect to some inertial frame vanishes identically. The fluid motion 

^ \ with respect to the inertial frame is then more or less a circular translation, whereas 

in a frame which follows the orbital motion (designed hereafter as the co- orbiting 
frame) ^ it is a counter-rotation. For a more extensive description of irrotational 
Darwin-Riemann configurations, we report to Eriguchi's review'^). 

The present article focuses on the general relativistic treatment of irrotational 
binary configurations. These configurations can be seen as generalizations of the 
irrotational Darwin-Riemann ellipsoids in the following directions: 

1. the fluid is no longer assumed to be incompressible; 

2. the gravitational potential of the companion is no longer truncated to the second 
order, but totally considered; 

3. the Newtonian treatment is replaced by a general relativistic one. 
The motivation for such a study is twofold: 



*' Contrary to MacLaurin ellipsoids, which are exact solutions for rotating incompressible fluids 
in Newtonian gravity, Darwin ellipsoids are approximate solutions, because of the second order 
truncation in the expansion of the gravitational potential of the companion. 
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1. Investigate the stability with respect to gravitational collapse: by 

means of numerical computations, Wilson, Mathews and Marronetti^^' have 
found that, due to relativistic effects, binary neutron stars may individually 
collapse to black hole prior to merger. This rather surprising result is now 
thought to be due to an error in some analytical formula implemented in the 
numerical code . Consequently this motivation is now rather weak. 

2. Provide realistic initial conditions for binary coalescence: it has been 
shown that the gravitational-radiation driven evolution of a neutron star binary 
system is too rapid for the viscous forces to synchronize the spin of each neutron 
star with the orbit ^^'^^ as they do for ordinary stellar binaries. Rather, the 
viscosity is negligible and the fluid velocity circulation (with respect to some 
inertial frame) is conserved in these systems. Provided that the initial spins 
are not in the millisecond regime, this means that close binary configurations 
are better approximated by zero vorticity (i.e. irrotationat) states than by 
synchronized states. These irrotational configurations constitute realistic initial 
conditions for fully hydrodynamical computations of the merging phase, as 
performed by different groups ' -"^^^ . 

The plan of this article is as follows. Having presented the general formalism 

to treat relativistically irrotational binary systems in Sect. 2, we specialize it to the 
case where the spatial 3-metric is assumed to be conformally flat in Sect. 3, and 
exhibit the full system of partial differential equations to be solved. Some analytical 
solutions are presented in Sect. 4. We then discuss numerical techniques to solve the 
problem in Sect. 5, and present the numerical results obtained by various groups in 
Sect. 6. The paper ends with a discussion about the innermost stable circular orbit 
in these numerical solutions (Sect. 7). 

§2. Formalism for relativistic irrotational binaries 

Generalizing the Newtonian formulation presented in Ref. 12) we have proposed 
a relativistic formulation for quasiequilibrium irrotational binaries The method 
is based on one aspect of irrotational motion, namely the counter-rotation (as mea- 
sured in the co-orbiting frame) of the fluid with respect to the orbital motion. This 
formulation has been slightly corrected by Asada-*^^) who has shown that in order 
to lead unambiguously to a counter-rotating state, the iterative procedure presented 
in Ref. 13) must be initiated with a vanishing velocity field with respect to the 
co-orbiting observer. Asada^^^ has also shown that the relativistic definition of 
counter-rotation implies that the flow is irrotational, i.e. that the vorticity 2-form 
vanishes identically. 

Subsequently, Teukolsky and Shibata gave independently two formulations 
based on the very definition of irrotationality, which implies that the specific enthalpy 
times the fluid 4- velocity is the gradient of some scalar fleld -"^^^ {potential flow) . 

The formulations presented by us^^^ (as amended by Asada^^^), Teukolsky 
and Shibata ^^-^ are equivalent; however the one given by Teukolsky and by Shibata 
greatly simplifles the problem. Consequently we used it in the following discussion. 

The general relativistic treatment of irrotational binary systems is based on 
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two physically well justified assumptions: (i) quasiequilibrium state (which means 
a steady state in the co-orbiting frame), and (ii) irrotational flow. Let us examine 
successively these two assumptions and their relativistic (geometrical) translation. 

2.1. Quasiequilibrium assumption 

When the separation between the centres of the two neutron stars is about 50 km 
(in harmonic coordinates) the time variation of the orbital period, Porbi computed 
at the 2nd Post-Newtonian (PN) order by means of the formulas established by 
Blanchet et al. is about 2%. The evolution at this stage can thus be still considered 
as a sequence of equilibrium configurations. Moreover the orbits are expected to 
be circular (vanishing eccentricity), as a consequence of the gravitational radiation 
reaction . In terms of the spacetime geometry, we translate these assumptions by 
demanding that there exists a Killing vector field I which is expressible as 

l = k + f2m, (2-1) 

where J7 is a constant, to be identified with the orbital angular velocity with respect 
to a distant inertial observer, and k and m are two vector fields with the following 
properties, k is timclikc at least far from the binary and is normalized so that far 
from the star it coincides with the 4-velocity of inertial observers, m is spacelike, 
has closed orbits, is zero on a two dimensional timelike surface, called the rotation 
axis, and is normalized so that V(m • m) • V(m • m)/(4m • m) tends to 1 on the 
rotation axis [this latter condition ensures that the parameter ip associated with m 
along its trajectories by m = d/dip has the standard 27r periodicity]. Let us call I 
the helicoidal Killing vector. We assume that iis a symmetry generator not only for 
the spacetime metric g but also for all the matter fields. In particular, I is tangent 
to the world tubes representing the surface of each star, hence its qualification of 
helicoidal (cf. Figure 1). 

The approximation suggested above amounts to neglect outgoing gravitational 
radiation. For non-axisymmctric systems — as binaries are — imposing I as an 
exact Killing vector leads to a spacetime which is not asymptotically flat ^'^^ . Thus, 
in solving for the gravitational field equations, a certain approximation has to be 
devised in order to avoid the divergence of some metric coefficients at infinity. For 
instance such an approximation could be the Wilson & Mathews scheme ^^-^ that 
amounts to solve only for the Hamiltonian and momentum constraint equations (cf. 
Sect. 3). This approximation has been used in all the relativistic studies to date and 
is consistent with the existence of the helicoidal Killing vector field (2T). Note also 
that since the gravitational radiation reaction shows up only at the 2.5-PN order, 
the helicoidal symmetry is exact up to the 2-PN order. 

Following the standard 3+1 formalism, we introduce a spacetime foliation by 
a family of spacelike hypersurfaces St such that at spatial infinity, the vector k 
introduced in Eq. (2-1) is normal to X!t and the ADM 3-momentum in Ut vanishes 
(i.e. the time t is the proper time of an asymptotic inertial observer at rest with 
respect to the binary system). Asymptotically, k = d/dt and m = d/dip, where if 
is the azimuthal coordinate associated with the above asymptotic inertial observer, 
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Fig. 1. Spacetime foliation Et, helicoidal Killing vector I and its trajectories = const, which are 
the worldlines of the co-orbiting observer (4- velocity: v). Also shown are the rotating-coordinate 
shift vector B and the unit future-directed vector n, normal to the spacelike hypersurface Et. 



SO that Eq. (2-1) can be re- written as 

I = 



d ^ d 

dt dif 



(2-2) 



One can then introduce the shift vector B of co-orbiting coordinates by means 
of the orthogonal decomposition of I with respect to the St fohation (cf. Fig. 1): 

l = Nn-B, (2-3) 

n being the unit future directed vector normal to X!t and N the lapse function. 

2.2. Irrotational flow 

We consider a perfect fluid, which constitutes an excellent approximation for 
neutron star matter. The matter stress energy tensor is then 

T = {e + p)u®u + pg , (2-4) 

e being the fluid proper energy density, p the fluid pressure and u the fluid 4- 
velocity. A zero-temperature equation of state (EOS) is a very good approximation 
for neutron star matter. For such an EOS, the first law of thermodynamics gives 
rise to the following identity (Gibbs-Duhem relation): 



Vp 



1 



e + p h 
where h is the fluid specific enthalpy: 

e+p 



V/i , 



h 



m-Qn 



(2-5) 



(2-6) 
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n being the fluid baryon number density and m-Q the mean baryon mass: me = 
1.66 X 10~^^ kg. Note that for our zero-temperature EOS, me h is equal to the fluid 

chemical potential. 

By means of the identity (2-5), it is straightforward to show that the classical 
momentum-energy conservation equation V • T = is equivalent to the set of two 
equations 22),23). 

u-(VAw) = 0, (2-7) 
V • (nw) = . (2-8) 
In Eq. (2-7), w is the co-momentum 1-form 

w:=hu (2-9) 

and V Aw denotes the exterior derivative of w, i.e. the vorticity 2-form^^). In terms 
of components, one has 

(V A w)a/3 = VaWp - VpWa ■ (2-10) 

The advantage of writing the equation of motion in the form (2-7) rather than in the 
traditional form V • T = is that one can see immediately that a flow of the form 

w = V^, (2-11) 

where ^ is some scalar, is a solution of the equation of motion. Such a flow is called 
a potential flow. Indeed, Eq. (2-11) implies the vanishing of the vorticity 2-form: 

VAmj = 0, (2-12) 

so that the equation of motion (2-7) is trivially satisfied. Equation (2-12) is the 
relativistic definition of an irrotational flow^^^. 

2.3. First integral of motion 

The above two assumptions, namely (i) Z is a symmetry generator and (ii) the 
flow is irrotational, yield to the following flrst integral 

/tZ • tt = const. (2-13) 

This was flrst pointed out by Carter . The demonstration is straightforward if one 
applies the classical Cartan's identity to the 1-form w and the vector fleld I: 

£lw = l-iV Aw) + V{l-w) , (2-14) 

where £i denotes the Lie derivative along the vector field I. Hypothesis (i) implies 
that = and hypothesis (ii) that V A lu = 0, so that Eq. (2-14) reduces to 

V(Z-iy) = 0, (2-15) 

from which the first integral (2-13) follows. 
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2.4. Remark about Bernoulli's theorem 

Let us mention a point which seems to have been missed by various authors: the 
existence of the first integral of motion (2-13) is not merely the relativistic general- 
ization of Bernoulli's theorem. This latter follows solely from the existence of the 
symmetry generator I and can be established as follows. Inserting £j^w = into 
Cartan's identity (2-14) yields 

I - {V ^w) + V{l■w) = Q . (2-16) 

Performing the scalar product by u leads to 

^w)■u + u■V{l■w) = Q . (2-17) 

The first term in the left-hand side vanishes by virtue of the equation of motion 
(2-7), so that one is left with 

M-V(Z-w) = 0, (2-18) 

which means that the quantity I w = hi u is constant along each streamline. This 
is the Bernoulli theorem. The key point is that, in order for the constant to be 
uniform over the streamlines (i.e. to be a constant over spacetime), as in Eq. (2-13), 
some additional property of the flow must be required. One well known possibility 
is rigidity (i.e. u colinear to I) ^'^^ which would apply to synchronized binaries. The 
alternative property with which we are concerned here is irrotationality. 

2.5. Equation for the velocity potential ^ 

Since Eq. (2-7) is trivially satisfied by the potential flow (2-11), the only part of 
the momentum-energy equation V • T = which remains to be satisfied is Eq. (2-8) 
(baryon number conservation). Inserting Eq. (2-11) in Eq. (2-8) results in an equation 
for the velocity potential 

Ti /" rL\ 

-V -VW^ + V^ -V i-j =0 . (2-19) 

Following the 3+1 formalism, we introduce the 3-metric h induced by g into the 
Et hyper surf aces, and denote by D the associated covariant derivative. Taking into 
account the helicoidal symmetry, Eq. (2-19) becomes*^ 

hr ( h B^ \ 

uDiD'^ + D'nDi^ = -^B'DiU + n I L>*<f A In — - ^ A^n J + RnhE,, , 

(2-20) 

where K is the trace of the extrinsic curvature tensor of the Et hypersurfaces and 
is the Lorentz factor between the fluid observer and the Eulerian observer (observer 
whose 4- velocity is the unit normal n to Et): 

1 . \V2 



Tn := -n • w = ( 1 + j^Di^ L>V ) . (2-21) 



*^ Latin indices . . . run from 1 to 3. 
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§3. A simplifying assumption: the conformally flat 3-metric 

3.1. Presentation and justifications 

The problem of finding irrotational binary configurations in quasiequilibrium 
amounts to solve the elliptic equation (2-20) for ^ and the Einstein equations for the 
components of the metric tensor. As a first step, a simplifying assumption can be 
introduced, in order to reduce the computational task, namely to take the 3-metric 
induced in the hypersurfaces Et to be conformally flat: 

h = A^f, (3-1) 

where A is some scalar field and / is a flat 3-mctric. Note that this assumption, first 
introduced by Wilson & Mathews , is physically less justified than the assumptions 
of quasiequilibrium and irrotational flow discussed above. However, some possible 
justifications of (3-1) are 

1. it is exact for spherically symmetric conflgurations; 

2. it is very accurate for axisymmetric rotating neutron stars 

3. the 1-PN metric fits it; 

4. the 2.5-PN metric deviates from it by only 2 % for two 1.4 Mq neutron stars as 
close as 30 km (in harmonic coordinates) . 

3.2. Partial differential equations for the metric 
Assuming (3-1), the full spacetime metric takes the form 

ds"^ = -{N^ - BiB')dt^ - 2Bidt dx' + A^fijdx' dx^ . (3-2) 

We have thus five metric functions to determine: the lapse N, the conformal factor 
A and the three components -B* of the shift vector. Let us introduce auxiliary metric 
quantities: the shift vector of non-rotating coordinates: 

N = B + f2—, (3-3) 



and the two logarithms 



(3 := In(^iV), (3-4) 
u:=lnN. (3-5) 



At the Newtonian limit (3 = and u coincides with the Newtonian gravitational 
potential. 

In the following, we choose maximal slicing coordinates, for which K = 0. 
The Killing equation Val/3 + V/jZq = give rise to a relation between the Ut 
extrinsic curvature tensor K and the shift vector N: 

where V stands for the covariant derivative associated with the flat 3-metric /. 
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The trace of the spatial part of the Einstein equation, combined with the Hamil- 
tonian constraint equation, result in the following two equations 

Al3 = ATrA'^S+-A'^KijK'^ --(Viu'Ti^ + ViPTp) , (3-7) 

Au = 4ttA^{E + S) + A^KijK'^ - VivV'p , (3-8) 

where A := V*Vi is the Laplacian operator associated with the flat metric /, E and 
S are respectively the matter energy density and trace of the stress tensor, both as 
measured by the Eulerian observer: 

E:=n-T ■n = rl{e + p)-p, (3-9) 

S :=h-T = ^p+{E+p)UiU\ (3-10) 
being the fluid 3-velocity as measured by the Eulerian observer: 

U:=^hu. (3-11) 

-1 n 

For the potential flow (2-11), U' is related to ^ by 

By the means of Eq. (3-6), the momentum constraint equation yields 
AN' + V* (VjNA = -16TrNA'^{E + p)W + 2N A^ K'^V 0^ - 4u) . (3-13) 

The equations to be solved to get the metric coeflficients are the elliptic equa- 
tions (3-7), (3-8) and (3-13). Note that they represent only five of the ten Einstein 

equations. The remaining five Einstein equations are not used in this procedure. 
Moreover, some of these remaining equations are certainly violated, refiecting the 
fact that the conformally flat 3-metric (3-1) is an approximation to the exact metric 
generated by a binary system. 

At the Newtonian limit, Eqs. (3-7) and (3-13) reduce to = 0. There remains 
only Eq. (3-8), which gives the usual Poisson equation for the gravitational potential 
v. 

3.3. Equation for the matter 

Taking the logarithm of the first integral of motion (2-13) and using the metric 
(3-2) yields 

iV 



+ + InTn + In 1 + A'^fij — U^ ) = const , (3-14) 



where H is the logarithm of the specific enthalpy: 

H:=lnh. (3-15) 
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At the Newtonian limit, H coincides with the specific non-relativistic (i.e. which 
does not include the rest mass energy) enthalpy. 

Taking into account Eqs. (2-21) and (3-3), the Newtonian limit of the first inte- 
gral (3-14) is 

H + V + ^iy^f -{nxr)-V^ = const. (3-16) 

(Recall that u reduces to the Newtonian gravitational potential). We recognize the 
classical expression [compare e.g. with Eq. (12) of Ref. 15) or Eq. (11) of Ref. 27)]. 

For a zero-temperature EOS, H can be considered as a function of the baryon 
density n solely, so that one can introduce the thermodynamical coefficient 

< - ^ . (3-17) 

dmn 

The gradient of n which appears in Eq. (2-20) can be then replaced by a gradient of 
H so that, using the metric (3-1), one obtains 

CHAIr + VHVi^ = -A^hr^^ViH + CH ^^Vi{H -(3)- A^^hV^^^'^ . 

(3-18) 

Note that this is not a Poisson-type equation for because the coefficient H in front 
of the Laplacian operator vanishes at the surface of the star. Numerically speaking, 
this means that this equation must be dealt by a specific technique and not by the 
direct use of some Poisson solver. 

At the Newtonian limit Eq. (3T8) reduces to 

(HA^ + VH -V^ = {n X r) -VH . (3-19) 

Here again, we recognize the classical expression [compare e.g. with Eq. (13) of 
Ref. 15)]. 



§4. Analytical approach 



The above equations cannot be solved analytically, unless additional simplifying 
assumptions are introduced. Two such assumptions are that 

1. the fluid is incompressible; 

2. the 1-PN approximation to general relativity is used. 

These assumptions have recently allowed Taniguchi^^) to find analytical solutions 
to the relativistic irrotational Darwin-Riemann problem. Note that at the 1-PN 

level, the figures are no longer ellipsoids, even for an incompressible fluid. This 
is explicitely taken into account in Taniguchi's procedure ^^-^ , which improves over 
a previous 1-PN study by Lombardi, Rasio & Shapiro where the figures where 
assumed to be ellipsoidal. 

The main findings of Taniguchi's study arc that the orbital separation at 
the innermost stable circular orbit (ISCO) *^ is lower than in the Newtonian case, 

*^ Taniguchi defines the ISCO as the location of the energy minimum along a constant baryon 
number sequence of decreasing separation, the true ISCO being certainly close to this point. 
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whereas the orbital angular at the ISCO is roughly the same than in the Newtonian 
case. 

§5. Numerical approaches 

Recently three groups have obtained numerical solutions of the partial differen- 
tial equations presented in Sect. 3, by the mean of different numerical techniques. 
In chronological order, these groups are 

• our group ^'^^'^^^ which has employed a multi-domain spectral method with 
spherical coordinates; 

• Marronetti, Mathews & Wilson > who have employed single-domain finite 
differences with Cartesian coordinates; 

• Uryu & Eriguchi , who have employed multi-domain finite differences with 
spherical coordinates. 

In this Section, we discuss only the numerical technique developed in our group, 
whereas in Sect. 6, we will present the results obtained by the three groups. 

Roche ellipsoid 



aj/a^ = 0.75 aj/a^ = 0.68 




5 10 15 20 25 30 35 



Fig. 2. Logarithm of the relative global error of the numerical solution with respect to the number of 
degrees of freedom in 6 for a Roche ellipsoid for an equal mass binary system and f2^/{iTGp) = 
0.1147 (the numbers of degrees of freedom in the other directions arc Nr — 2Ng — 1 and 
A^,^ = Ne — 1) Also shown is the error in the verification of the virial theorem. 

5.1. Brief description of the multi-dom,ain, spectral m,ethod, 

The numerical procedure to solve the PDE system presented in Sect. 3 is based 
on the multi-domain spectral method developed in Ref. 31). We simply recall here 
some basic features of the method: 

• Spherical-type coordinates 0' , ip') centered on each star are used: this ensures 
a much better description of the stars than by means of Cartesian coordinates. 

• These spherical-type coordinates are surface-fitted coordinates: i.e. the surface 
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of each star lies at a constant value of the coordinate ^ thanks to a mapping 
{^,6', if') i—f {r,9,ip) [see Ref. 31) for details about this mapping]. This en- 
sures that the spectral method applied in each domain is free from any Gibbs 
phenomenon (spurious oscillations generated by discontinuities). 

• The outermost domain extends up to spatial infinity, thanks to the mapping 
1/r = (1 — S)/{2Rq). This enables us to put exact boundary conditions on 
the elliptic equations (3-7), (3-8) and (3-13) for the metric coefficients: spatial 
infinity is the only location where the metric is known in advance (Minkowski 
metric) . 

• Thanks to the use of a spectral method in each domain, the numerical error 
is evanescent, i.e. it decreases exponentially with the number of coefficients (or 
equivalently grid points) used in the spectral expansions, as shown in Fig. 2. 

The PDE system to be solved being non-linear, we use an iterative procedure. 
The iteration is stopped when the relative difference in the enthalpy field between 
two successive steps goes below a certain threshold, typically 10~^. An illustrative 
solution is shown in Fig. 3. 




Fig. 3. Velocity field with respect to the co-orbiting frame in a irrotational binary system (only the 
part of the stars above the orbital plane is represented). 

The numerical code is written in the LORENE language which is a C-|— |- 
based language for numerical relativity. A typical run make use of A^r = 33, A'g = 21, 
and = 20 coefficients (= number of collocation points, which may be seen as 
number of grid points) in each of the domains. 8 domains are used : 3 for each star 
and 2 domains centered on the intersection between the rotation axis and the orbital 
plane. The corresponding memory requirement is 155 MB. A computation involves 
~ 250 steps, which takes 9 h 30 min on one CPU of a SGI Origin200 computer 
(MIPS RIOOOO processor at 180 MHz). Note that due to the rather small memory 
requirement, runs can be performed in parallel on a multi-processor platform. This 
especially useful to compute sequences of configurations. 
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5.2. Tests passed by the code 

In the Newtonian and incompressible limit, the analytical solution constituted 
by a Roche ellipsoid is recovered with a relative accuracy of ~ 10~^, as shown in 
Fig. 2. For compressible and irrotational Newtonian binaries, no analytical solution 
is available, but the virial theorem can be used to get an estimation of the numer- 
ical error: we found that the virial theorem is satisfied with a relative accuracy of 
10~^. Preliminary comparisons with the irrotational Newtonian configurations re- 
cently computed by Uryu & Eriguchi^^^' ^'''^ reveal a good agreement. Regarding the 
relativistic case, we have checked our procedure of resolution of the gravitational 
field equations by comparison with the results of Baumgarte et al. which deal 
with corotating binaries [our code can compute corotating configurations by setting 
to zero the velocity field of the fluid with respect to the co-orbiting observer]. We 
have performed the comparison with the configuration za = 0.20 in Table V of 
Ref. 38). We used the same equation of state (EOS) (polytrope with 7 = 2), same 
value of the separation rc and same value of the maximum density parameter q™^^. 
We found a relative discrepancy of 1.1% on 17, 1.4% on Mq, 1.1% on M, 2.3% on J, 
0.8% on ZA, 0.4% on and 0.07% on rs (using the notations of Ref. 38)). 

§6. Numerical results 

6.1. Equation of state and compactification ratio 

As a simplified model for nuclear matter, let us consider a polytropic EOS with 
an adiabatic index 7 = 2: 

p = K{mBny , e = uiBn + p/{'y — 1) , (6-1) 

with K = 1.8 X 10~^ Jm^kg"^. This EOS is the same as that used by Mathews, 
Marronetti and Wilson (Sect. IV A of Ref 39)). 

The mass - central density curve of static configurations in isolation constructed 
upon this EOS is represented in Fig. 4. The three points on this curve corresponds 
to three configurations studied by the various groups mentioned in the beginning of 
Sect. 5: 

• The configuration of baryon mass Mb = 1.625 Mq and compactification ratio 
M/R = 0.14 is that considered in the dynamical study of Mathews, Marronetti 
and Wilson and in the quasiequilibrium studies of our group of Mar- 
ronetti et al. and of Uryu &: Eriguchi ^"^^ . 

• The configuration of baryon mass Mb = 1.85 M© and compactification ratio 
M/R = 0.17 has been studied by our group and Uryu &: Eriguchi "^'^^ . 

• The configuration of baryon mass Mb = 1.95 Mq and compactification ratio 
M/R = 0.19 has been studied by Marronetti et al.^^)'-'^^) *). 

*^ Marronetti et al. use a different value for the EOS constant ft: their baryon mass 

Mb = 1.55 Mq must be rescaled to our value of k in order to get Mb = 1.95 Mq. Apart from this 

scaling, this is the same configuration, i.e. it has the same compactification ratio M/R = 0.19 and 
its relative distance with respect to the maximum mass configuration, as shown in Fig. 4, is the 
same. 
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Spherical models in isolation 

Polytrope y=2 k=0.0332 p„„^"'c'' 

2.2 
2 
1.8 

s" 

— 1.6 
(/) 
(/) 

(0 

1.4 

1.2 
1 

2 4 6 8 10 12 14 16 18 
Central energy density e^ [ p^j^ ] 

Fig. 4. Mass as a function of the central energy density for static isolated neutron stars constructed 
with the EOS (6-1). The heavy dots are configurations considered by our group Mar- 
ronetti at al. and Uryu & Eriguchi^*' (see text) (p„uc := 1-66 x 10^'' kgm"^). 

6.2. Irrotational sequence with M/R = 0.14 

In this section, we give some details about the irrotational sequence Mb = 
1.625 Mq presented in Rcf. 30). This sequence starts at the coordinate separation 
d = 110 km (orbital frequency / = 82 Hz), where the two stars are almost spherical, 
and ends at d = 41 km (/ = 332 Hz), where a cusp appears on the surface of the 
stars, which means that the stars start to break. The shape of the surface at this 
last point is shown in Fig. 5. 

The velocity field with respect to the co-orbiting observer, as defined by Eq. (52) 
of Ref. 13), is shown in Fig. 6. Note that this field is tangent to the surface of the 
star, as it must be. 

The lapse function N is represented in Fig. 7. The coordinate system (x, y, z) 
is centered on the intersection between the rotation axis and the orbital plane. The 
X axis joins the two stellar centers, and z = is the orbital plane. The value of N 
at the center of each star is Nc = 0.64, whereas the central value of the conformal 
factor is = 2.20. The shift vector of non-rotating coordinates, N, is shown in 
Fig. 8. Its maximum value is 0.10 c. 

The variation of the central density along the Mb = 1.625 Mq sequence is shown 
in Fig. 9. We have also computed a corotating (i.e. synchronized) sequence for com- 
parison (dashed line in Fig. 9). In the corotating case, the central density decreases 
quite substantially as the two stars approach each other. This is in agreement with 
the results of Baumgarte et al. . In the irrotational case (solid line in Fig. 9) , the 
central density remains rather constant (with a slight increase, below 0.1%) before 
decreasing. This contrasts with the result of the dynamical calculations by Mathews 
et al. which showed a central density increase of 14% for the same compactifica- 
tion ratio M/R = 0.14. But, as stated in Sect. 1, this last result has revealed to be 
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Fig. 5. Shape of irrotational binary neutron stars of baryon mass A4b = 1.625 A/q, when the 
coordinate separation between their centers (density ma^xima) is 41 km. Only one half of each 
star is represented (the part which is above the orbital plane). The drawing is that of the 
numerical grid, which coincides with the surface of the star, thanks to the use of surface- fitted 
spherical coordinates. 



Velocity / co— orbiting observer 




-10 

X [km] 



Fig. 6. Velocity field with respect to the co-orbiting observer, for the configuration shown in Fig. 5. 
The plane of the figure is the orbital plane. The heavy line denotes the surface of the star. The 
companion is located at x = +41 km. 

spurious due to some error in a formula employed in the code^^. 

It is worth to note that, for the same compactification ratio M/R = 0.14, the 
computations by the other groups show a central density variation in accordance 
with the one presented above: it is at most 1% both in the results of Marronetti et 



Darwin- Riemann problems in GR 



15 




lapse N 




-50 50 

x [km] 



Fig. 7. Isocontour of the lapse function A'^ for the configuration shown in Fig. 5. The plots are cross 
section in the x = 0, j/ = and z = Q planes (note that the x coordinate is shifted by 20.5 km 
with respect to that of Fig. 6). The dot-dashed line denotes the boundary between the inner 
numerical grid and the outer compactified one (which extends to spatial infinity), for the grid 
system centered on the intersection between the rotation axis and the orbital plane. 



al. and Uryu & Eriguchi '^'^^ . 

We can thus conclude that no tendency to individual gravitational collapse is 
found when the orbit shrinks. 

6.3. Irrotational sequence with M/R = 0.17 

In order to investigate the dependency of the above result on the compactness 
of the stars, we have computed an irrotational sequence with a baryon mass Mb = 
1.85 Mq, which corresponds to a compactification ratio M/R = 0.17 for stars at 
infinite separation (second heavy dot in Fig. 4) . The result is compared with 
that of M/R = 0.14 in Fig. 10. A very small density increase (at most 0.3%) is 
observed before the decrease. For the same compactification ratio M/R = 0.17, 
Uryu & Eriguchi report a slightly higher increase (1.5 %). 

This density increase remains within the expected error (~ 2%, cf. Sect. 3.1) 
induced by the conformally flat approximation for the 3-metric, so that it cannot be 
asserted that this effect would remain in a complete calculation. 
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shift N°i dans le plan orbital 




X [km] 

Fig. 8. Shift vector of non-rotating coordinates in the orbital plane, for the configuration shown in 
Figs. 5-7. 



Model 



Existence of an ISCO References 



Newtonian corotating 

Newtonian irrotational 
GR corotating 
GR irrotational 



ISCO <S4. 7 > 2 42) 
ISCO « 7 > 2.4 37) 

ISCO <^ 7 > 5/3 40), 42) 
ISCO exists for 7 = oo 28) 



Table I. Known results about the existence of an ISCO for extended fluid bodies, in terms of the 
adiabatic index 7. 



6.4. Irrotational sequence with M/R = 0.19 

Marronetti, Mathews and Wilson have computed quasiequihbrium irro- 

tational configurations with a compactification ratio M/R = 0.19 (third heavy dot 
in Fig. 4). They found a central density increase, as the orbit shrinks, of the or- 
der of 1.5%. This still remains within the error introduced by the conformally fiat 
approximation. 



§7. Innermost stable circular orbit 



An important parameter for the detection of a gravitational wave signal from 
coalescing binaries is the location of the innermost stable circular orbit (ISCO), if 
any. In Table I, we recall what is known about the existence of an ISCO for extended 
fluid bodies. The case of two point masses is discussed in details in Ref. 
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Fig. 9. Relative variation of the central energy density Cc with respect to its value at infinite sepa- 
ration ej;"' as a function of the coordinate separation d (or of the orbital frequency J?/(27r)) for 
constant baryon mass Mb = 1.625 Mq sequences The solid (resp. dashed) line corresponds 
to a irrotational (resp. corotating) sequence of coalescing neutron star binaries. 



0.006 



Irrotational sequences 

EOS: Y=2 K=0.0332 p„„^"'c' 



0.004 
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Separation d [ l<m ] 



Fig. 10. Relative variation of the central energy density ec with respect to its value at infinite 
separation ejf' as a function of the coordinate separation d for constant baryon mass sequences 
with Mb = 1.625 M0 (solid line, same as in Fig. 9) and Mb = 1.85 M© (dashed line). 



For Newtonian binaries, it has been shown that the ISCO is located at a 
minimum of the total energy, as well as total angular momentum, along a sequence 
at constant baryon number and constant circulation (irrotational sequences are such 
sequences). The instability found in this way is dynamical. For corotating sequences, 
it is secular instead This turning point method also holds for locating ISCO 
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in relativistic corotating binaries For relativistic irrotational configurations, no 
rigorous theorem has been proven yet about the locahzation of the ISCO by a turning 
point method. All what can be said is that no turning point is present in the 
irrotational sequences considered in Sect. 6: Fig. 11 shows the variation as the orbit 
shrinks of the ADM mass of the spatial hypersurface t = const (which is a measure 
of the total energy, or equivalently of the binding energy, of the system), as well as of 
the total angular momentum, for the evolutionary sequences with M/R = 0.14 and 
M/R = 0.17. Clearly, both the ADM mass and the angular momentum decreases 
monoticaly, without showing any turning point. The same result has been found by 
Marronetti et al. [Fig. 2 of Ref. 33) , which also shows the good agreement between 
Marronetti et al. results and ours 26)] and Uryu & Eriguchi [Fig. 5 of Ref. 34)]. 



Irrotational sequence U^=^ .625 




40 50 60 70 SO 90 100 110 120 130 140 
Separation d [km] 

Irrotational sequence Mb=1.85 M^, 

EOS: 1-^ K.0.0332 
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Irrotational sequence Mg=1 .625 M^^ 

EOS: 7-2 rc.0.0332 
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Irrotational sequence Mb=1 .85 M^, 

EOS: 7-2 1C.0.0332 




40 50 60 70 80 90 100 110 120 
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Fig. 11. Variation of (half of) the ADM mass (left) and of the total angular momentum (right) of 
the binary system with respect to the coordinate distance d, along the evolutionary sequences 
Mb = 1.625 Mq, M/R = Q.U (top) and Mb = 1.85 M©, M/7? = 0.17 (bottom) [from Ref.^^)] 
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